A new approach for efficient simulation of Coulomb interactions in ionic fluids 
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We propose a simplified version of local molecular field (LMF) theory to treat Coulomb inter- 
actions in simulations of ionic fluids. LMF theory relies on splitting the Coulomb potential into 
a short-ranged part that combines with other short-ranged core interactions and is simulated ex- 
plicitly. The averaged effects of the remaining long-ranged part are taken into account through a 
self-consistently determined effective external fleld. The theory contains an adjustable length pa- 
rameter a that specifies the cut-off distance for the short-ranged interaction. This can be chosen 
to minimize the errors resulting from the mean-field treatment of the complementary long-ranged 
part. Here we suggest that in many cases an accurate approximation to the effective field can be 
obtained directly from the equilibrium charge density given by the Debye theory of screening, thus 
eliminating the need for a self-consistent treatment. In the limit cr — > 0, this assumption reduces to 
the classical Debye approximation. We examine the numerical performance of this approximation 
for a simple model of a symmetric ionic mixture. Our results for thermodynamic and structural 
properties of uniform ionic mixtures agree well with similar results of Ewald simulations of the full 
ionic system. In addition we have used the simplified theory in a grand-canonical simulation of a 
nonuniform ionic mixture where an ion has been fixed at the origin. Simulations using short-ranged 
truncations of the Coulomb interactions alone do not satisfy the exact condition of complete screen- 
ing of the fixed ion, but this condition is recovered when the effective field is taken into account. We 
argue that this simplified approach can also be used in the simulations of more complex nonuniform 
systems. 



I. INTRODUCTION 

The long-ranged nature of the Coulomb inter- 
action often causes problems in computer simula- 
tionSfi Although in some cases direct truncation of 
Coulomb interactions, e.g., by reaction field meth- 
ods^ or shifted force truncation^ can give accu- 
rate results^ such methods have suffered from sig- 
nificant errors in many other physically relevant 
casesj^i^ii Truncation schemes tend to work best in 
dense uniform systems, where there is considerable 
cancellation of the long-ranged electrostatic forces, 
and they often perform poorly in inhomogcneous 
systems i^i^i^i^ At the same time, the Ewald sum 
methodiSiii — which does not truncate Coulomb 
interactions and accurately accounts for all images 
generated by periodic boundary conditions — is 
generally complex and computationally demanding 
when applied to inhomogcneous systemsji^ although 
considerable simplification is possible in reduced ge- 
ometries such as a two-dimensional slab.— i^^'^'i^'^ 
Finally, the strict periodicity of the Ewald sum 
method has been known to introduce artifacts in 
some biologically relevant studies.^ 

In this paper, we examine the performance of 
an alternative treatment of inhomogcneous Coulomb 
systems, local molecular field (LMF) theory^i^ when 
applied to a simple model for an ionic solution. This 



theory accounts for the averaged effects of slowly- 
varying long-ranged components of the intermolec- 
ular interactions by using a self-consistently deter- 
mined effective field. LMF theory provides a general 
framework for treating both uniform and nonuni- 
form systems. We show for the ionic system con- 
sidered in this paper that analytic results from the 
Debye theory of screening can be used to greatly 
simplify the theory and find very good agreement 
with results of conventional simulations using Ewald 
sums. 

The derivation of LMF theory and its application 
to general Coulombic systems are given in detail else- 
where^i^i^i^ii^iS^ The ideas behind the theory are 
best understood by considering a nonuniform one 
component system with long-ranged intermolecular 
interactions w{r) in an external field (t>{r), which can 
represent the interactions with fixed objects such as 
walls or solutes. LMF theory relates structural and 
thermodynamic properties of the original system to 
those of a "mimic system" with short-ranged inter- 
actions uo{r) in a renormalized effective field 4>r(t) 
that accounts for the averaged effects of the remain- 
ing long-ranged component ui(r) = w{r) — uo{r) of 
the intermolecular interactions. A key idea in LMF 
theory is that ui{r) should be properly chosen to be 
slowly varying so that the averaging can give accu- 
rate results. 

For such a ui{r), the effective field is determined 
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in principle by the condition that the nonuniform 
singlet density in the mimic system (denoted by the 
subscript R) equals that in the original system: 



[<Pr]) = [0]). 
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An explicit equation for (?!)fl,(r) can be derived by 
subtracting the first equations of the exact Yvon- 
Born-Green hierarchy that relate the gradient of the 
singlet density to forces in the full and mimic sys- 
tems. As argued in Refs. HiHOlHlillll , when ui 
is chosen to be slowly varying over the range of pair 
correlations between neighboring molecules, the ef- 
fective field 4>H{r) is accurately given by the self- 
consistent LMF equation. 

M^) = J dv'pniv'- [(A«])wi(|r'-r|)+C, 

where C is a constant of integration. Equation ([T]) 
then relates structure in the original and mimic sys- 
tems; thermodynamic properties can be similarly re- 
lated by integration over the structure. Finally, uni- 
form systems are treated in LMF theory by choosing 
t/)(r) to be the field arising from a fixed fluid particle, 
i.e., by taking 4>{y) = w(r). This field will induce a 
nonuniform singlet density that can be directly re- 
lated to the pair distribution function in the uniform 
fluid34 

LMF theory has many ideas in common with ear- 
lier methods, but they arc implemented in new ways 
that avoid most of the limitations of those methods. 
As can be seen from Eq. LMF theory uses a 
mean-field average of the long-ranged interactions, 
similarly in spirit to the random-phase approxima- 
tion of density functional theory^^ or to the Debye 
theory of screening for ionic systems.— But, crucially, 
the average in LMF theory is taken only over partic- 
ular slowly varying components ui , chosen precisely 
so that their averaged effects can be accurately de- 
scribed by an effective field. 

The general idea of separating intermolecular in- 
teractions into short-ranged and long-ranged parts 
has also proved useful in many different contexts, 
particularly for uniform systems where long-ranged 
forces largely cancel and thus can be treated as a 
weak perturbation or by standard integral equation 
closures.— i^i^'^'^i^ For historical reasons the in- 
termolecular potential w was often split so that the 
properties of the reference system with the inter- 
molecular potential uq would be well-known, e.g., a 
hard sphere fluid or an ideal gas. However, these 
choices do not generally guarantee that the remain- 
ing interactions ui can in fact be accurately treated 
by perturbation or diagrammatic methods. 

With present day computers, complex short- 
ranged systems can be simulated efficiently and so 



the nature of the reference system is no longer an 
overriding issue. Therefore the split of the inter- 
molecular potential w in LMF theory can be op- 
timized to minimize the errors associated with the 
approximate mean-field treatment of its long-ranged 
part ui. The mimic system is then defined as the 
special short-ranged reference system resulting from 
this optimal split of w, and we can use the simula- 
tions to accurately determine its properties. 

LMF theory thus corrects the two major short- 
comings of the classical Debye theory of ionic sys- 
tems, namely, its inaccurate mean-field averaging of 
the entire Coulomb potential and the highly approx- 
imate Boltzmann form of the density response to 
the effective field. But as discussed below the (lin- 
earized) Debye theory satisfies the exact Stillinger- 
Lovett moment conditions^ and therefore correctly 
describes the asymptotic behavior of the charge cor- 
relation function. These features of the Debye the- 
ory can be exploited to greatly simplify the deter- 
mination of the effective field in LMF theory while 
still giving accurate results, as we now show. 



II. APPLICATION OF LMF THEORY TO A 
SYMMETRIC IONIC FLUID MODEL 

We consider a uniform mixture of N positive and 
N negative ions with charges +q and —q. The molec- 
ular cores are described by the repulsive part of the 
LJ potential u^"^ {r)i^ The total intermolecular po- 
tential between ions with charges qi and qj (= ±g) 
is taken to be 



(3) 



where r is the distance between the ion centers. 

The strength of the Coulomb interactions can be 
characterized by the ratio F = Ig/d of the Bjer- 
rum length Ib = q^/ksT (the distance where the 
Coulomb energy between two positive ions equals 
ksT) to the "ion diameter" d = a^'^ , where ct^'^ is 
the length parameter in the repulsive LJ potential, 
r > 1 characterizes the "strong coupling" regime. 

Previous studies^i^ i^^'^^ have shown that it is ad- 
vantageous to divide the Coulomb interaction v{r) = 
1/r into short-ranged and long-ranged parts in the 
following way (see Fig. [T|): 



1 crfc(r/cr) crf(r/cr) 



VQ{r,a) + vi{r,a), (4) 



r r r 

where cr is a length scale at our disposal. The idea 
behind this separation is best understood by consid- 
ering the Fourier transform of the long-ranged com- 
ponent Vi{r,a) = CT:i{r/(T)/r, 



vi{k,(j) 



fc2 



exp 
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FIG. 1: Separation of the Coulomb potential v{r) = 
(solid) into a short-ranged part vo{r,a) — r~^eric{r / a) 
(dash) and a long-ranged part vi{r,a) — r^^eri{r / a) 
(dash-dot). The cut-off length a = 1.5d shown here is 
also used in Fig. 2(a) as discussed below in Sec. III. 



The function vi{k^a) differs significantly from zero 
only at small wave- vectors fca < 2. As a result, 
vi (r, a) remains finite as r — > and is slowly varying 
for r < a, while still decaying asymptotically as 1/r. 
This makes it much more suitable for a mean-field 
averaging than the full Coulomb potential v{r) used 
in the classical Debye theory. 

This generates a separation of the intermolecu- 
lar potentials Wij{r) = UQ^ij(r) + (r) into short- 
ranged and long-ranged components respectively, 
where 



and 



UQ,tjir) = UQ-\r) + q.qjVoir, a) (6) 



UiAjir) = qiqjVi{r,a). (7) 



As mentioned above, LMF theory for a uniform 
system focuses on the density response to a fixed 
fiuid particle. For the symmetric ionic system con- 
sidered here, we can assume without loss of general- 
ity that a positive ion is fixed at the origin. This 
yields a single particle field (f>j{r) ~ w-^-j{r) act- 
ing on an ion with charge qj. Using Eqs. © and 
([7]), 0j(r) naturally separates into a short-ranged 
core part uo^+j^r) and the long-ranged remainder 
ui,+j{r). 

The fixed ion induces a nonuniform singlet density 
Pj{r) = pg^j{r), proportional to the radial distri- 
bution function g^j{r) in the uniform fluidj^ Here 
p ~ N/V is the number density of positive or neg- 
ative ions. Similarly, the induced charge density 
p^(r) = qp_|_(r) — qp^{r) satisfies 

p'(r) = qp [g++ir) - ,9+_(r)] ee qph^r). (8) 



The electrostatic energy per ion in the uniform ionic 
fluid can be written exactly in terms of the charge 
density as 



^ = 1 1 dv'p'^{v')v{r'). 



(9) 



LMF theory models the original uniform system 
by a nonuniform mimic system comprised of "sol- 
vent" mimic ions with short-ranged intermolecu- 
lar interactions UQ^ij{r) in an effective field (jiRji^), 
which we can picture as arising from a modified "so- 
lute" ion fixed at the origin. According to Eq. ([T]), 
when LMF theory is accurate, the solute-induced 
densities in the mimic system should equal the anal- 
ogous densities in the original system. Using Eqs. 
([6]) and ([7]), the LMF equation for (j)iij{r) can be 
written as [cf. Eq. ©] 

q, f dr'[q5{r') + p%{r')]v,{\r-v'\,<j),{10) 



where p%{r) — qpR^+{r) — qpR-{r) is the induced 
charge density in the nonuniform mimic system. 

The form of the effective field in Eq. ([T0|) clearly 
depends on a. There are two criteria that help deter- 
mine an optimal choice of a. For the LMF method 
to be quantitatively valid, a should be large enough 
that vi (r, a) remains slowly varying on the scale of 
short-ranged pair correlations. At the same time, it 
is desirable to keep a small in order to reduce simula- 
tion times of the mimic system. Thus a is generally 
chosen near its minimal accurate value fTmin, which 
is state dependent and of the order of a characteristic 
neighbor spacing . ^'^'^'^ 

In principle, Eq. (fTO|) has to be solved self- 
consistently since the effective field (fiuj (r) depends 
on the charge density /3^(r) in the presence of the 
field itself. During an iterative solution of Eq. (|10p . 
the density induced by a given field can be accurately 
determined from the simulation of the nonuniform 
mimic system. This procedure has been successfully 
carried out for models of iona^ or water— confined 
between charged and uncharged hard walls. Once 
the self-consistent charge density p|j (r) has been de- 
termined, it can be used in Eq. ([9]) to calculate the 
electrostatic energy of the original ionic system. 

However, self-consistent simulations of the 
nonuniform mimic system can be computationally 
demanding and may not always be required to ob- 
tain accurate results. In the following we introduce 
a simplified version of LMF theory that avoids the 
need for full self-consistency. We show that both 
the local structure and the electrostatic energy 
of the uniform ionic system can be accurately 
reproduced by combining results of straightforward 
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(a)r = 5, p<f = 0.3816 
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(b)r = 5, p(f = 0.0012 

FIG. 2: The ion-ion distribution functions gij{r) in two 
uniform symmetric ionic mixtures with F = 5, pd'^ — 
and F = 5, pd^ = 0.0012 fCb)! The open 



0.3816 (a) 



5, pd-" = 0.0012 [( 

symbols are the ion pair distribution functions in the 
full Coulomb systems, obtained using Ewald simulations. 
The curves are the same functions in the mimic systems 
with[(a)]o-/d = 1.5 (dash) and 1 (dash-dot), [(b)]a-/d = 15 
(dash) and 10 (dash-dot). The g+-{r) functions in (a) 
are vertically displaced by 0.5. The insets show a zoom 
of the first peak of g+-{r). 



simulations of the uniform mimic system with 
(/)flj(r) = along with analytic results from the 
Dcbye theory of screening for corrections arising 
from the long-ranged forces. Similar ideas have 
also proved useful in simulations of more complex 
systems of charged polymers, as will be described 
elsewhere 1^ 



III. LOCAL STRUCTURE FROM 
SHORT-RANGED SIMULATIONS 

A. Strong coupling approximation 

At high density in the uniform ionic system we 
expect considerable cancellation of the forces from 
the slowly- varying part vi (r, a) of the Coulomb in- 
teraction when a is chosen properly. In the sim- 
plest "strong-coupling approximation" (SCA) to the 
full LMF theoryj^iiSii^ we ignore all effects of vi {r, a) 
on the fluid structure, i.e., neglect the integral in 
Eq. pI7|) . Note that the strong short-ranged part 
Wo (r, a) of the Coulomb interaction as well as the 
LJ core is still taken into account in the SCA. This 
part of the interaction would be expected to domi- 
nate local structural arrangements at lower densities 
as well, provided that a is chosen large enough. 

In the SCA the self-consistent field 0ijj(r) in 
Eq. (jlOp is approximated by (/)o.j(r) = uo.+j(r), the 
known field of a solvent mimic ion. In this case, 
the induced charge density p^ir) can be determined 
directly and more efficiently from the radial distribu- 
tion functions gQ_+j{r) in the uniform mimic system 
[cf. Eq. dSD], 

PoW = QP [ffo,++W - 5o,+- W] = 9P^oW- (11) 

Thus the SCA can be viewed as a particularly use- 
ful direct truncation scheme^ whose accuracy can 
be justified in certain limits and corrected, if neces- 
sary, by the full LMF theory. Note that we use the 
subscript to refer both to pair correlation functions 
in the uniform mimic system with — and to 
the equivalent singlet densities in the nonuniform 
system when (j^nj is approximated using the SCA 
by (ko-j- The subscript R refers to the nonuniform 
mimic system in the presence of the full renormal- 
ized field given by Eq. pUj) . 

Figure [2] gives representative results of Langevin 
dynamics simulations using the Ewald sum method 
for the present ionic fluid model at states with fixed 
temperature ksT/e^'^ = 1 and charge q chosen so 
that the system is at moderately strong coupling 
with F = 5. These are compared to simulations 
of the uniform mimic system, as described by the 
SCA, where the cut-off radius for the short-ranged 
potential vo{r^a) is 2.5a. 



Figure 2(a) shows results for a high density state 
with pd^ = 0.3816, where there is substantial cancel- 
lation of attractive forces. We find that for a/d = 1.5 
there is excellent agreement between the distribution 
functions in the full Coulomb and mimic systems 
over the range of r shown. Higher values of cr give 
equally good results, but smaller a values cause no- 
ticeable errors, as illustrated in the inset for a/d = 1, 
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so (Tmin is about 1.5o? for this state. This good agree- 
ment is consistent with previous work on LMF the- 
or3*^i2ii2^ and with earher findings that models with 
truncated Coulomb interactions can often provide a 
good description of structural features in dense uni- 
form systems.—!^ 

Figure 2(b) [ makes the same comparison for a very 
low density state with pcfi = 0.0012. Despite this 
low density the coupling is strong enough that the 
Debye theory alone would give poor results, and this 
state presents a major challenge to theory. We find 
that a much larger value of a/d — 15 is needed to 
achieve comparably good results (the small devia- 
tions for cr/d = 10 are shown in the inset). Large a 
is needed since there is essentially no force cancella- 
tion at low densities and the characteristic neighbor 
distances are large, but this also makes simulations 
of the mimic system much more costly in this regime. 

However, despite the excellent agreement between 
the distribution functions in the original and uni- 
form mimic systems in the range of r shown, there 
are fundamental differences in the asymptotic be- 
havior of these functions as r ^ oo. This is most 
easily seen from the small- fc behavior of the charge 
structure factor S"^(fc), which is simply related to the 
Fourier transform of the dimensionless charge corre- 
lation function /i'(7-) defined in Eq. ([5]), 



^^(fc) = l + p/i'?(fc). 



(12) 



As A: — *■ 0, the charge structure factor S'^{k) of any 
ionic system exhibits the same universal behavior. 



where 



Xd = (STr^Bp)-!/^ 



(13) 



(14) 



is the Debye screening length. The exact form in 
Eq. is independent of any details of the short- 
ranged core interactions Mo,ij(r) and is a conse- 
quence of the Stillinger-Lovett moment conditions 
In contrast, the analogous function 



(15) 



for the uniform short-ranged mimic system, where 
/iQ(r) is defined in Eq. (jlip . will remain finite as 
/c — > 0, with the coefficient of k^ depending on the 
details of the intermolecular interactions. 

In Fig. [3] we compare SQ{k) and 5"'(fc) for the same 
ionic mixtures whose gij{r) are shown in Fig. [21 For 
o" ^ cTmin, the functions S'g(fc) closely follow S"'(fc) 
everywhere except for small fc, where the differences 
described above can be seen. If cr < CTmin, noticeable 
discrepancies between S'o(fc) and S"'(fc) appear also 
at larger fc, as illustrated by the dash-dot curve in 
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(a)r = 5, pd^ = 0.3816 
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FIG. 3: The charge structure factors S'^{k) for the same 
systems as in Fig. (2] The insets show a zoom of the small- 
k region where a discrepancy is seen in the structure 
factors of the full Coulomb and mimic systems. The solid 
curves show the results of application of the perturbation 
equation (|27p to the dashed curves in Figs. [3tiO] and [3(b)] 



B. Complete screening and the Debye theory 

This different behavior at small wave-vectors im- 
plies that the uniform mimic system will not exactly 
satisfy the basic "complete screening condition" that 
true ionic fluids obey. Complete screening (equiva- 
lent to the Stillinger-Lovett zeroth moment condi- 
tion) requires that the exact p'^ induced by a fixed 
positive ion in a grand ensemble will satisfy 



drp«(r) 



(16) 



Fig. 3(a) 



However, the results above for the simpler SCA im- 
ply that Eq. (fTB)) will not hold if is approximated 
by the p^ given by a grand canonical simulation of 
the mimic system with a fixed positive "solvent" 
mimic ion at the origin. 
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To verify these conclusions, we have performed 
grand-canonical ensemble simulations of the nonuni- 
form mimic system in the field (j)oj(r) = uo,+j(r) 
for a state with F = 5, pd^ = 0.0012 and a = 
lOd. We find that the mean number of counteri- 
ons (N) = 1202.19 in the simulation box exceeds the 
mean number of coions by AA^ = 0.68±0.10. This is 
definitely smaller than the mean difference A7V = 1 
that should hold in the case of complete screening. 
We chose a very large simulation box with L = lOOd 
such that this value of AA^ is independent of L and 
carried out a long 50ns simulation run to obtain the 
reported statistical convergence of AA^. 

In addition, we ran a similar simulation of the 
nonuniform mimic system in which the long-ranged 
part Vi (r, a) of the Coulomb potential has been 
taken into account through the LMF equation (fTU|) . 
One can show that the density induced by a self- 
consistent solution of this equation will exactly sat- 
isfy the complete screening condition.— However, in- 
stead of solving Eq. PH)) self-consistcntly, we have re- 
placed the charge density profile p'j^{r) in this equa- 
tion by the screening profile of a point charge given 
by the linearized Debyc-Hiickcl theory: 



■ cxp 



A, 



(17) 



At first glance this may seem to be a very crude 
approximation, since the Debye profile is generally 
accurate only when both T and p are very small. 
Otherwise p%i{r) and p%{r) will differ considerably 
at small r. However, the Debye profile has the cor- 
rect asymptotic behavior since it satisfies the exact 
Stillinger-Lovett zeroth and second moment condi- 
tions.— Moreover, when integrated over the slowly 
varying part of the Coulomb potential vi (r, a) as in 
Eq. (dni), most of the short-ranged features of this 
profile on the scale r < cr become irrelevant, so the 
resulting estimate for (j)ji j{r) can still be accurate. 

With this approximation, the integration in 
Eq. (fTO|l can be carried out exactly and we obtain 
an explicit expression for the effective field 4>Rj{r), 



^R,j{r) « uo,+j{r) + _ cxp ( ^ 



exp 



exp 



D 



r \ ^ /(J r 

erfc 

V2 cr 



A 



D 



r \ ^ f(7 r 



(18) 



The simulation of the nonuniform mimic system, 
where (/)_Rj(r) is given by Eq. (fTS]) . yields AA^ = 
1.09 ± 0.10. Thus our simple estimate for the ef- 
fective field using the Debye theory can reproduce 
the complete screening behavior seen in the full sys- 
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pcP = 


0.0012 


pd^ = 


0.3816 




cr = 10 


(7 = 15 


cr = 1 


cr = 1.5 




-0.5946(4) 


-0.6837(3) 


-0.4912(1) 


-1.3122(2) 


Debye 


-0.8487(4) 


-0.8622(3) 


-3.2068(1) 


-3.1597(2) 


RPA 


-0.8453(5) 


-0.8612(3) 


-3.0136(2) 


-3.1335(2) 


Debye-M 


-0.8678(4) 


-0.8701(3) 


-3.2979(1) 


-3.1906(2) 




-0.8708(4) 


-3.1880(3) 



TABLE I: The electrostatic energy U for the same sys- 
tems as in Figs. [2]and|3l The error in the last significant 
figure is indicated in parentheses. Approximations for 
the total energy labeled Debye, RPA, and Debye-M are 
discussed in Eqs. p3)l . p8)l . and ([29j respectively. 



tem or from a complete self-consistent solution of 
the LMF equation. 



IV. ELECTROSTATIC ENERGY 

Thermodynamic properties of ionic systems also 
require careful attention to contributions from the 
long-ranged parts of the Coulomb interactions. 
Again we find that analytic results from the Debye 
theory can provide simple but accurate corrections 
to results from the uniform mimic system. In anal- 
ogy with Eq. ([2]), the "electrostatic energy" of the 
uniform mimic system is given by 



= dv'pl{v')v,{r',a) 



(19) 



As illustrated in TablclH t/o differs considerably from 
the full Coulomb energy U , determined by the Ewald 
sum method, even for a > (Tmin when local structural 
properties of the original and mimic systems closely 
resemble each other. 

To find the needed correction, we note that Eq. 
^ can be exactly rewritten as 



2N 



2 (27r) 



J dk[S'>{k)-l]Mk,>j) 
(20) 



+ ^ / dr'[h'^{r')-hl{r')]vo{r',a), 



f3U 

m 



where h'^{r), h1^{r). S'^{k) and Uq arc defined in 
Eqs. JHl), (HI]), and respectively. 

We expect that the value of the last integral in 
Eq. (|20p is very small since with proper choice of a 
h'^{r) and h1{r) are very similar over the entire range 
of r where (r, a) differs significantly from zero (see 
Fig. [2]). Hence the energy difference AU = U — Uq 
can be accurately estimated as 



(3AU Ib 1 



2N 



2 (27r)= 



■ j dkvi{k,cr) [S''{k) - 1] 



(21) 
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Note that the LMF integral in Eq. (fTU)) for r = 
equals AC//iV. 

The function vi{k, ct), given by Eq. ([5]), is a rapidly 
decaying function of k for ka > 2. Thus only the 
small-fc behavior of S'^{k) is significant in Eq. ((2T|) . 
Similar to our discussion of complete screening in 
Sec. IIII B[ this suggests that we can accurately use 
the Debye approximation 



Slik) 



fc2 + • 



(22) 



for the charge structure factor S'^{k) in Eq. ((2T|) . 

Equation is exact at small enough k since 
it satisfies both Stillingcr-Lovctt moment conditions 
and, unlike Eq. (|13p . correctly reduces to unity at 
large k. Furthermore it becomes an exact result for 
all k in the limit of very small T and p. Substituting 
Eq. ^ in Eq. (HJ), we find 



/?AC/ (3Ud 



2N 



2N 



h 



a 



(23) 



where Ud is the well known result for the Coulomb 
energy in the Debye approximation, 



and 



2N 



fi{y) = exp 



2\i 



y- 1 erfc 



(24) 



(25) 



We expect accurate results from Eq. (|23p only 
when a is properly chosen to be greater than a state- 
dependent minimum value dmin- For strong coupling 
states with F > 1, we note that (Tmin ^ Xd- Using 
the asymptotic expansion of erfc(y/2) in Eq. (pS]) . 
Eq. then reduces to the strong couphng energy 
correction 



I3MJ_ 
2N 



Ib 

./ttct 



D 



(26) 



derived in Ref . [23l. where S''{k) was approximated by 
the second moment term in Eq. But Eq. 

also correctly reduces to the exact Debye energy 
in the limit of very weak coupling and low density 
where (Tmin ^ and provides a more generally useful 
expression. 

The resulting energy estimates C/o + At/ from 
Eq. ((23l) , labeled "Debye" , are given in Table |T] and 
are plotted for more states over a wide range of a in 
Fig. m We see that the deviations of the Debye cor- 
rected energy from the Ewald energies U e\y are re- 
duced by approximately an order of magnitude from 
the uncorrected results given in Table HI Moreover, 
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(b)pd3 



: 0.02984 



FIG. 4: The Debye corrected energies Uo + AU from 
Eq. (|23|l for a wide range of a relative to the Ewald 
energies Uew of the full Coulomb systems. The data 
sets, shown from left to right, are for |(a)| F = 5 and 
pd^ = 0.0012, 0.005969, 0.02984, 0.2387, 0.3816, and for 
Wi\pd ^ = 0.02984 and F = 0.2, 1, 5, 10. The two stars 
in |(a)| show the location of (Tmin as estimated from the 
structure shown in Fig. [5] for a low density state with 
pd'^ — 0.0012 and (Jmin ~ 15, and a high density state 
with pd"^ — 0.3816, and (Tmin ~ 1.5. A spline curve is 
fitted through each set of data to guide the eye. 



the accuracy of the correction AU sharply increases 
with larger values of the ratio a/Xo- This allows us 
to determine appropriate values for (Tmin and obtain 
energy estimates accurate to within 1% when a is 
chosen large enough. 

In the limit a s- 0, Uo becomes vanishingly small 
and Eq. (|23p reduces to the classical Debye approx- 
imation. Since the Debye theory generally overesti- 
mates the absolute value of U, all the curves in Fig. 2] 
turn up sharply at small enough a < (Tmin- 

These results emphasize that LMF theory is an in- 
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herently approximate approach and accurate results 
can be expected only when a is chosen greater than 
some state-dependent minimum value CTmin- In the 
present simplified treatment, errors for a wide range 
of a have been determined from comparison with 
the results of Ewald simulations. These fairly basic 
simulation studies can be used to devise an estimate 
of (Tniin for an ionic solution of given concentration 
and ionic strength. 

Although the energy correction in Eq. ([23|) de- 
pends only on the ratio ct/Ad, the numerical per- 
formance of this correction depends significantly on 
another length scale, the the characteristic neigh- 
bor distance between a pair of ions, estimated here 
by r* = (2p)-i/3. As can be deduced from Fig. [H 
the accuracy of AJ7 decreases sharply when a < r* . 
When a/r* > 1, the accuracy of AJ7 is generally ac- 
ceptable and is higher for larger values of cr/ Ad. As 
a rule of thumb for systems with moderately strong 
coupling, wc suggest that if an ionic solution is sim- 
ulated as part of a more complex system, the two 
conditions a > 2r* and cr > SA^ should be fulfilled 
simultaneously to yield an accurate description of 
electrostatic interactions. 



INTEGRAL EQUATION METHODS FOR 
A UNIFORM IONIC FLUID 



Integral equation methods that treat the long- 
ranged part of intcrmolccular interactions as a weak 
perturbatio n^^i^^ have long been used to study ionic 
fluids and could serve as an alternative approach to 
LMF theory for uniform systems. These equations 
are usually derived by summing certain classes of 
diagrams, where individual diagrams represent dif- 
ferent terms in the perturbation series for structural 
and thermodynamic quantities. Although it is dif- 
ficult to develop a physical intuition for what the 
errors will be in a given application, these meth- 
ods are expected to work best when the effects of 
the long-ranged perturbations on the structure of 
a short-ranged reference system are in some sense 
small. 

Since the local structure is well described by the 
uniform mimic system, it seems likely that this could 
serve as a particularly useful reference system for 
perturbation integral equation methods. This idea 
was in fact suggested long ago by Ceperley and 
Chester)^ although they did not discuss the possi- 
bility of applying it outside the framework of integral 
equations. Here we use one of the earliest perturba- 
tion approaches, the RPA-like method of Ref. H^, 
to correct results for the uniform mimic system. In 
this approach the charge structure factor S"'(fc) of a 
uniform ionic mixture can be approximately written 



as 



S'i{k) = 



Sl{k) 



l + 2lBpvi{k,a)Sl{ky 



(27) 



where S^^k) is the charge structure factor of the uni- 
form mimic system and vi{k, a) is given in Eq. ([5]). 

The resulting functions S'^{k) arc shown in Fig. [3] 
where they are compared with S^i^k) obtained from 
short-ranged simulations. In contrast with S'o(fc), 
the functions S'^{k) satisfy both Stillinger-Lovett 
moment conditions. For cr > CTmin they are virtually 
indistinguishable from the results of Ewald simula- 
tions. 

The perturbation method of Ref. |2^ also yields an 
expression for the energy correction AC/, 



(3AU Ib 1 



2N 



2 (27r)3 



1 



■In 



Sl{k) 



vi{k,a) 



21bp"' S^{k) 

(28) 

where S'^{k) is given by Eq. p7|) . The energy esti- 
mates resulting from Eq. (|28p . labeled as "RPA" are 
given in Table [H 

Remarkably, although the calculation based on 
the Debye profile is much less involved, it gives en- 
ergy estimates that are equally accurate. We find, 
however, that yet another energy correction gives 
even more accurate results, 



f3AU I3Ud PUoD 



2N 



2N 



2N 



(29) 



Here Uqd is the "electrostatic energy" of the uni- 
form mimic system in the Debye limit, obtained by 
summing for "Coulomb cores" voir, a) the same ring 
diagrams that lead to the conventional Debye ex- 
pression Ud when using the full Coulomb interaction 
v{r). This gives 



2N 2N ■'^ V 



(30) 



where 



f3{y) 



1 — exp — 



p + l_expl 



dk. (31) 



The idea behind Eq. (|29p is that with a proper 
choice of ct, the energy correction A[/ should be in- 
dependent of most details of the short-ranged in- 
teractions. Most errors in the Debye treatment of 
the short-ranged part of the Coulomb interactions 
are canceled by subtraction of the two terms in Eq. 
((29)) . The results in Table U for this "Debye-Mimic" 
(Debye-M) approximation give best agreement with 
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the Ewald energies Uew of the full Coulomb sys- 
tems. This seems to indicate that approximations 
of the RPA type work best if perturbation terms, 
similar to At/, contain no information on the short- 
ranged core structure of the mimic system. 

VI. CONCLUDING REMARKS 

In this paper we have used a simplified version 
of local molecular field (LMF) theoryi^ to calculate 
the structural and thermodynamic properties of a 
symmetric ionic mixture. LMF theory has already 
been applied successfully to the description of uni- 
form and nonuniform ionic systemSf^'^'2^ where the 
LMF was either approximated using only the SCA or 
was determined by a full self-consistent calculation. 
We have shown for the ionic mixture considered in 
this paper that we can go beyond the SCA but avoid 
the necessity of finding a self-consistent solution if we 
replace the mimic system's charge density with its 
Debye analogue in the LMF equation (fTO|l . The re- 
sulting effective field, given in Eq. (fT5|) . is sufficiently 
accurate to reproduce the exact complete screen- 
ing condition in grand-canonical ensemble simula- 
tions. Furthermore, the energies of uniform ionic 
mixtures, obtained under this assumption, agree well 
with those calculated using diagrammatic techniques 
and the Ewald sum method. 

When applied to complex inhomogencous sys- 
tems, this assumption will significantly speed up 
the simulations of underlying mimic systems since 
it reduces the simulations in an a priori unknown 
self-consistent field to the simulations in a simpler 
and well-defined external field. In addition, the De- 
bye charge density is found analytically from the 



Debye-Hiickel equation for an infinite system. For 
unbounded Coulomb systems this should result in 
more accurate values of the effective field in Eq. ^ 
than those obtained by numerical integration using 
simulation results in a finite volume. 

These ideas are being actively applied in our stud- 
ies of polyelectrolytes in salt solutionsi^ Since the 
polymer dynamics are slow in comparison with the 
dynamics of small ions, it is valid to assume that the 
distribution of salt ions around a polyelectrolyte is 
always in local equilibrium. In this case, the total 
effective field of a polyelectrolyte is simply a sum of 
the effective fields of all its monomer charges, given 
by Eq. p^ . This essentially amounts to replac- 
ing the full Coulomb potential of a monomer charge 
with a screened Coulomb potential, defined as a sum 
of r~^erfc(r/i7) and the last term in Eq. p8|). We 
note that the idea of introducing effective Debye- 
Hiickel, or Yukawa, interactions between polyelec- 
trolyte charges has been widely used to account for 
the screening by small ionsi^i^ The advantage of 
our approach, however, is that it invokes only the 
long-ranged features of the Debye screening profile, 
while still explicitly accounting for the strong elec- 
trostatic core interactions between charges at small 
distances in the simulations. This approach should 
therefore remain accurate for strongly interacting 
and dense polyelectrolyte systems, where conven- 
tional Debye-Hiickel interactions are a very crude 
approximation. 
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